import numpy as np
import pandas as pd

def get_ngc1068_diffuse_equivalent():

    # energies in TeV
    e_range=np.logspace(np.log10(1.5), np.log10(15), 50)

    norm_unit = 1e-14  # convert to GeV-1 already
    norm_at_1tev = 5.0 * norm_unit
    gamma = 3.2

    flux = norm_at_1tev * e_range**(-gamma)
    e_in_gev = e_range * 1e3
    esq_phi_conv = (e_in_gev)**2 / (4 * np.pi)

    # also get uncertainties from contour
    gammas = ngc_1068_contour_68CL[::2]
    norms = ngc_1068_contour_68CL[1::2] * norm_unit
    unc_lo, unc_up = get_spl_butterfly(norms, gammas, energies=e_in_gev)

    return e_in_gev, flux * esq_phi_conv, unc_lo * esq_phi_conv, unc_up * esq_phi_conv


def get_dnn_fermi_pi0():
    # for relative uncertainties
    temp_bf = 21.8
    rel_unc_up = (temp_bf + 5.3) / temp_bf
    rel_unc_lo = (temp_bf - 4.9) / temp_bf

    # get bestfit flux line
    x, f = dnn_fermi_pi0_e, dnn_fermi_pi0_esq_phi / (4 * np.pi)
    f_lo = f * rel_unc_lo
    f_up = f * rel_unc_up
    return x, f, f_lo, f_up


dnn_fermi_pi0_all = np.array([899.228436403271, 5.869268556669518e-7,
    999.8625546729646, 5.480080589484773e-7,
    1106.1569532602514, 5.105948636722797e-7,
    1223.7514041579564, 4.755809117353856e-7,
    1353.84720474313, 4.432568284193927e-7,
    1497.773361128012, 4.1286057026992135e-7,
    1657.0001647492713, 3.844234380656769e-7,
    1833.1542122709984, 3.5829510756706e-7,
    2028.0350222387524, 3.337250843095814e-7,
    2243.6334182336213, 3.1094125559764983e-7,
    2482.1518663212037, 2.8961851431429605e-7,
    2746.026974554672, 2.6984589627573224e-7,
    3037.9543843777415, 2.5142317958938146e-7,
    3360.916308208025, 2.3418187977949193e-7,
    3718.2119945136496, 2.181228990366823e-7,
    4113.491431602015, 2.0323137603330075e-7,
    4550.792634424944, 1.8935651592198383e-7,
    5034.582895307215, 1.7637142653086657e-7,
    5569.804419999234, 1.6427678733458059e-7,
    6161.924815253234, 1.5306140569276962e-7,
    6816.99294368379, 1.4261171217654445e-7,
    7541.700716503276, 1.3287543230034699e-7,
    8343.451455381808, 1.2372319905959307e-7,
    9230.435521788912, 1.152764615873477e-7,
    10211.713986415685, 1.0744139954685172e-7,
    11297.311193410258, 1.0004101912029429e-7,
    12561.611681186776, 9.296841338425947e-8,
    13827.000892804803, 8.731567801079349e-8,
    15143.169915673861, 8.102381516554645e-8,
    16753.025321668247, 7.588685495355205e-8,
    18627.883659488736, 7.081916475965468e-8,
    20608.19553464137, 6.575372644855812e-8,
    22799.03240525542, 6.071922024246782e-8,
    25350.509238752813, 5.658940787511146e-8,
    28045.496785612206, 5.2853768764685844e-8,
    31026.985791252002, 4.922933097480988e-8,
    34325.433942194206, 4.5883332519544316e-8,
    37974.53685147017, 4.2736891155075344e-8,
    42011.57227938255, 3.981919016843207e-8,
    46477.78094803687, 3.7088596286824996e-8,
    51418.78783988818, 3.45452523948519e-8,
    56885.06828411976, 3.2176317857787045e-8,
    63251.16888518717, 2.9902432676208144e-8]
)

dnn_fermi_pi0_e = dnn_fermi_pi0_all[::2]
dnn_fermi_pi0_esq_phi = dnn_fermi_pi0_all[1::2]


ngc_1068_contour_68CL = np.array(
    [
        2.8671311243509523, 2.915749907735755,
        2.873839202502996, 3.2565708188251143,
        2.889880258953535, 2.6651447928685936,
        2.885505425376115, 3.4755595024149155,
        2.8986299261083746, 3.6856902167808574,
        2.913212704699774, 3.9041867909693337,
        2.9219623718546135, 2.6418814393513923,
        2.9292537611503127, 4.12268336515781,
        2.946753095459992, 4.352334419109662,
        2.9540444847556917, 2.674986980895101,
        2.9642524297696715, 4.572143285035006,
        2.98612659765677, 2.730461131589962,
        2.98175176407935, 4.790311786289264,
        3.000709376248169, 5.015510421848173,
        3.0182087105578486, 2.802040680873656,
        3.0211252662761283, 5.255941015067151,
        3.0415411563040875, 5.487935447120552,
        3.050290823458927, 2.8843571625499056,
        3.0619570463320462, 5.710087691147445,
        3.0823729363600054, 2.974726343520567,
        3.0823729363600054, 5.923803774008761,
        3.1042471042471043, 6.140331910591934,
        3.1144550492610836, 3.0749377125177393,
        3.127579549993343, 6.358090320678421,
        3.1465371621621623, 3.183201780809327,
        3.152370273598722, 6.575575336653062,
        3.1654947743309814, 3.2467286307986054,
        3.1800775529223806, 6.79385319555206,
        3.210701387964319, 7.006444456924632,
        3.2427835008653974, 7.2006039843566505,
        3.274865613766476, 7.356289504048684,
        3.3069477266675547, 7.471711527268642,
        3.339029839568633, 7.5468700540165194,
        3.347779506723472, 4.057367026436438,
        3.371111952469711, 7.5763966180960445,
        3.379861619624551, 4.229157944717304,
        3.4031940653707897, 7.552238520212798,
        3.4119437325256294, 4.427791193979557,
        3.4352761782718684, 7.472606271634687,
        3.4425675675675675, 4.638414017746827,
        3.4673582911729466, 7.318710240674745,
        3.4688165690320867, 4.849775005616083,
        3.490690736919185, 5.059097254536992,
        3.495065570496605, 7.10855715769965,
        3.509648349088004, 5.290623010970084,
        3.514023182665424, 6.892275075817139,
        3.5271476833976836, 5.546109808491519,
        3.5271476833976836, 6.669349517016734,
        3.534439072693383, 6.453313489834885,
        3.536772317268007, 5.803729080085363,
        3.538813906270803, 6.221529961810624,
        3.53832781365109, 5.995979819536483,
        3.256756756756757, 3.6416938110749157,
    ]
)

def get_spl_butterfly(norms_at_1_tev, gammas, energies=np.logspace(4, 7, 121)):


    def spl(e_in_tev, phi_at_1_tev, gamma):
        return phi_at_1_tev * (e_in_tev)**(-gamma)

    # use energy in TeV
    flux_range = spl(
        energies[:, np.newaxis]/1e3, norms_at_1_tev[np.newaxis],
        gammas[np.newaxis]
    )

    return flux_range.min(axis=1), flux_range.max(axis=1)